Dune workflow on the baron dataset

1 Load data

We download the data from segerstolpe et al, hosted [here]

suppressPackageStartupMessages({
  library(SingleCellExperiment)
  library(stringr)
  library(scRNAseq)
})
sce <- BaronPancreasData()
# We select the largest batch to avoid pre-processing issues
# sce <- sce[, sce$donor == "GSM2230759"]
filt <- rowSums(counts(sce) >= 2) >= 10
sce <- sce[filt, ]
sce

## class: SingleCellExperiment 
## dim: 12336 8569 
## metadata(0):
## assays(1): counts
## rownames(12336): A1CF A2M ... ZZZ3 pk
## rowData names(0):
## colnames(8569): human1_lib1.final_cell_0001 human1_lib1.final_cell_0002
##   ... human4_lib3.final_cell_0700 human4_lib3.final_cell_0701
## colData names(2): donor label
## reducedDimNames(0):
## altExpNames(0):

rm(filt)

2 Pre-processing

The dataset can be normalized using SCVI

suppressPackageStartupMessages({
  library(Seurat)
})
pre_process_time <- system.time({
  se <- CreateSeuratObject(counts = counts(sce),
                           min.cells = 0,
                           min.features = 0,
                           project = "de")
  se <- AddMetaData(se, as.data.frame(colData(sce)))
  se <- NormalizeData(se, verbose = FALSE)
  se <- FindVariableFeatures(se, selection.method = 'vst', nfeatures = 4000,
                             verbose = FALSE)
  se <- se[VariableFeatures(se), ]
  se <- ScaleData(object = se, vars.to.regress = c("nCount_RNA", "donor"))
  sce <- as.SingleCellExperiment(se)
})

suppressPackageStartupMessages({
  library(reticulate)
})
use_python("/usr/local/bin/python3")
scvi <- import("scvi")
anndata <- import("anndata")
np <- import("numpy")
sc <- import("scanpy")
scvi_time <- system.time({
  adata <- anndata$AnnData(X = as.sparse(t(counts(sce))),
                           obs = data.frame(cells = colnames(sce),
                                            batch = sce$donor))
  scvi$data$setup_anndata(adata, batch_key = "batch")
  model <- scvi$model$SCVI(adata)
  model$train(n_epochs = as.integer(60), n_epochs_kl_warmup = as.integer(30))
})

We can visualize the data using the labels from the original publication

library(scater)
reducedDim(sce, "scvi") <- model$get_latent_representation()
denoised <- t(model$get_normalized_expression(adata, library_size = 10e4))
dimnames(denoised) <- dimnames(counts(sce))
assay(sce, "denoised") <- log1p(denoised)
sce <- runTSNE(sce, dimred = "scvi")
plotTSNE(sce, colour_by = "label")

plotTSNE(sce, colour_by = "donor")

3 Creating inputs to Dune

3.1 sc3

suppressPackageStartupMessages(library(SC3))
sc3_time <- system.time({
  sce_sc3 <- sce
  logcounts(sce_sc3) <- assay(sce, "denoised")
  rowData(sce_sc3)$feature_symbol <- rownames(sce_sc3)
  counts(sce_sc3) <- as.matrix(counts(sce_sc3))
  logcounts(sce_sc3) <- as.matrix(logcounts(sce_sc3))
  sce_sc3 <- sc3_estimate_k(sce_sc3)
  K <- metadata(sce_sc3)$sc3$k_estimation
  # Note: with R >= 4.0, RStudio and Mac OS, this can fails. 
  # A workaround is running 
  # parallel:::setDefaultClusterOptions(setup_strategy = "sequential")
  sce_sc3 <- sc3(sce_sc3, ks = K, n_cores = NCORES, rand_seed = 786907,
                 svm_num_cells = round(.5 * ncol(sce)))
  sce_sc3 <- sc3_run_svm(sce_sc3, ks = K)    
  sce$SC3 <- colData(sce_sc3)[, paste0("sc3_", K, "_clusters")] %>% as.factor()
})

plotTSNE(sce, colour_by = "SC3")

3.2 Seurat

seurat_time <- system.time({
  se <- RunPCA(se, verbose = FALSE)
  se <- FindNeighbors(se, verbose = FALSE)
  se <- FindClusters(object = se, verbose = FALSE)
  sce$seurat <- Idents(se)
})
plotTSNE(sce, colour_by = "seurat")

3.3 Seurat with SCVI

seurat_scvi_time <- system.time({
  seu <- as.Seurat(x = sce, counts = "counts", data = "counts")
  seu <- FindNeighbors(seu, reduction = "scvi", verbose = FALSE)
  seu <- FindClusters(object = seu, verbose = FALSE)
  sce$seurat_scvi <- Idents(seu)
})
plotTSNE(sce, colour_by = "seurat_scvi")

4 Running Dune

library(Dune)
df <- colData(sce)[, c("SC3", "seurat", "seurat_scvi")] %>%
  as.matrix()
dune_time <- system.time(
  merger <- Dune(clusMat = df, BPPARAM = BPPARAM, metric = "NMI")
)

plotPrePost(merger)

NMItrend(merger)

plotNMIs(merger$initialMat)

plotNMIs(merger$currentMat)

5 Picking the final clustering result to use

5.1 Visual pick

sce$seurat_Final <- as.factor(merger$currentMat$seurat)
plotTSNE(sce, colour_by = "seurat_Final")

sce$seurat_scvi_Final <- as.factor(merger$currentMat$seurat_scvi)
plotTSNE(sce, colour_by = "seurat_scvi_Final")

sce$SC3_Final <- as.factor(merger$currentMat$SC3)
plotTSNE(sce, colour_by = "SC3_Final")

5.2 Picking based on sihouette

library(cluster)
dist_mat <- dist(as.matrix(reducedDim(sce, "scvi")))
sils <- lapply(merger$currentMat %>% as.data.frame, function(label){
  silhouette(label, dist = dist_mat)[,3] %>%  mean()
}) %>% unlist()
sils

##         SC3      seurat seurat_scvi 
## -0.03325991  0.18725471  0.31371226

6 Runtimes

times <- c(pre_process_time[1],
           scvi_time[1],
           sc3_time[1],
           seurat_time[1],
           seurat_scvi_time[1],
           dune_time[1])
names(times) <- c("Pre-Process",
                  "SCVI",
                  "Sc3",
                  "Seurat",
                  "Seurat after SCVI",
                  "Dune")
df <- data.frame(times = times,
                 Name = factor(names(times), levels = names(times)))
ggplot(df, aes(x = Name, y = times, fill = Name)) +
  geom_col() +
  theme_classic() +
  labs(x = "Step", y = "Time (second)") +
  scale_color_brewer(palette = "Dark2") +
  scale_y_log10()